< 

in 



C/2 



Hubbard model with bond-charge interaction on a triangular 
lattice: a renormalization group study 

Bibhas Bhattacharyya^ and Shreekantha Si^ 



p • ^ Department of Physics, Scottish Church College, 1 & 3 Urquhart Square, Calcutta 700 

006, India 



Institut fiir Theoretische Physik, Universitat zu Koln, Ziilpicher Strasse 77, D-50937, 

Koln, Germany 



Abstract 

We have studied the Hubbard model with bond-charge interaction on a triangular lattice for a 
half-filled band. At the point of particle-hole symmetry the model could be analyzed in detail 
in two opposite regimes of the parameter space. Using a real space renormalization group we 
calculate the ground state energy and the local moment over the whole parameter space. The 
' RG results obey the exact results in the respective limits. In the intermediate region of the 

Q . parameter space the RG results clearly show the effects of the non-bipartite geometry of the 

O i lattice as well as the absence of symmetry in the reversal of the sign of the hopping matrix 

element. 

^ _ PACS No.: 71.10.Fd, 64.60.Ak, 71.10.Hf 
^ ' I. Introduction 

Q I In spite of an extensive effort over the last few years the problem of electronic cor- 

^ ' relation in low dimensional systems remains yet to be clearly understood. Models of 
' interacting electrons are difficult to handle in one and two dimensions owing to strong 
On ■ fluctuations. While there are a few exact solutions in one dimension (1-D) [1] and in infi- 
nite dimensions [2] the situation is worse in two dimensions (2-D) . Standard techniques 
like mean field approximation or perturbative calculations are of very limited use in treat- 
ing the intermediate to strong correlation in low dimensional systems which are typically 
^ • dominated by strong fluctuations. Numerical simulations and exact diagonalizations are 
P ■ also limited by small cluster size because the dimension of the Hilbert space soon becomes 
too large to be handled as one goes from one to two dimensions. Therefore, it seems en- 
terprising to apply an approximate real space renormalization group (RG) known as the 
^ I block RG (BRG) [3,4] in this context. This works reasonably well for 1-D systems [5- 7] 
' and over the whole range of the coupling constants. This method employs a truncation 
of the Hilbert space, retaining only a few low lying states only, to bring out the essential 
ground state properties. Although the efficacy of the truncation procedure is much more 
satisfactorily handled in a recently developed RG scheme known as the density matrix 
RG (DMRG) [8], the latter is yet to be developed for a truly 2-D system. The BRG has 
already been applied to interacting electrons in 2-D [9] but only for bipartite lattices. It 
is interesting to see how does it work in case of a non-bipartite lattice like the triangular 
one. On the other hand the problem of interacting electrons on a triangular lattice has 
been addressed for a long time [10- 13] in view of a rich phase structure including the 
possibility of frustration [14]. Particularly, ^He on graphite was considered to be a good 
example of the Hubbard model [15] on a triangular lattice [10]. Recent works on the or- 
ganic superconductors like k-BEDT-TTF compounds also investigate the phase diagrams 
of interacting-electron models on anisotropic triangular lattices [16]. In the present work 



we study a generalized Hubbard model with bond-charge interaction on the 2-D trian- 
gular lattice at half- filling by using the BRG. In Sec. II we present the model and the 
renormalization scheme as suitably adapted to the present problem. Sec. Ill summarizes 
some of the exact solutions of this problem in limiting cases. In Sec. IV we focus on our 
results obtained from the RG and make comparisons with the exact solutions. Sec. V 
concludes the present work. 

II. The Model and the RG scheme 

The generalized Hubbard model with bond-charge interaction is defined as follows: 

(ij),o- i 

+UY^ni^nii, (1) 

i 

where cl^{cifj) creates (annihilates) a particle with spin o" in a Wannier orbital located 
at the site i; the corresponding number operator is rijo- = cj^Cio-. Here rij = rijo- -|- rij-o-. 
Sum over (ij) denotes contributions from distinct nearest neighbour pairs of sites on a 
triangular lattice, /i is the chemical potential. We focus our attention to the special case 
with X = t for which the present model becomes particle-hole (p-h) symmetric even on 
the triangular lattice. In fact it has been shown that for some realistic systems X ~ t 
[17]. Further the choice X = t results in some simplification in implementing the RG 
and forces /i = f//2 for half-filling. Moreover, at this special value of X we have some 
exact results for comparison. One should remember, however, that there is no t — >■ —t 
symmetry due to the non-bipartite nature of the lattice. 

Now, ior X — t and for a half-filled band {{rii) = 1) one can put = U/2 to rewrite 
the Hamiltonian (2) as follows: 

H ^ -tYl cLcj>(l - rii-a - rij-a) + t^IZ - nil) (t; - na) 

{ij),a i ^ ^ 

+DY.h + ^N, (2) 

i 

where, the constant D = will account for the renormalization of the vacuum energy. 

Here we have added a constant term NU/2 to the Hamiltonian given by (1) to compensate 
for the chemical potential term subtracted thereof. We just keep aside this constant term 
for the purpose of renormalizing the parameters in the Hamiltonian and then add up 
the same to the ground state energy at the end. The Hamiltonian has several conserved 
quantities: total number of particles u, total spin S and the ^-component Sz of the total 
spin besides the p-h symmetry pointed out earlier. We divide the 2-D lattice with N sites 
into triangular blocks of three sites each (Fig. 1). The block-Hamiltonian is then 
diagonalized for the three-site block. Since we are interested in the ground state properties 
of the system we truncate the Hilbert space for the three-site blocks and retain only four 
low-lying states of the block-Hamiltonian governed by the symmetries mentioned above 
together with the point group symmetry of the block. It is easy to see that the block- 
Hamiltonian is block-diagonalizable in terms of the good quantum numbers u, S and Sz- 
The states with {iy^2,S^Sz^ 0}, ^ 3,S ^ ^ 1/2}, {u^3,S ^ -S^ = 1/2} 
and {u — A, S — Sz — 0} are considered for this purpose. Of these the states in the first 



and the fourth group are connected by the p-h symmetry while those in the second and 
third by spin-reversal symmetry. However, the point group symmetry in the 2-D lattice 
imposes further restriction on the choice of the states [9]. 

To resolve this point we consider the symmetry group C^^ [18] of the basic triangular 
block. We take our eigenstates of the block-Hamiltonian to also be the eigenstates of 
i?(C3„), the matrix representations of the group elements of C^^. While choosing the 
four states to be retained, we simply take one from each of the four groups so that all 
of them lie in the same representation of C^^ such that the contribution to the ground 
state energy be the minimum. This leads to the following RG equations which relate the 
renormalized parameters (primed quantities) to the original parameters (unprimed ones) 
in the Hamiltonian : 

U' = 2{E2-E,), 
D' = 3D + {E2 + E3)/2, 
t' = 2Re[X*{X-2X')]t , (3) 

where E2 and E-^ are the lowest eigenvalues in the subspaces {u — 2, S — Sz — 0} and 
{ly — 3, S = Sz = 1/2} respectively; and, 

X={iy = 2, S = Sz = Oh^liy = 3, S = Sz = l) 

X'^{u^2, S^Sz^ 0\cb^nbi\iy ^3, S^Sz^^) . (4) 

A* is the complex conjugate of A and Re denotes the real part. Here the subscript b refers 
to the site index of a boundary site of the block. For a triangular block, however, this 
could be any of the three sites. The factor of 2 appears in the renormahzation of the 
hopping due to the two connecting paths between the two neighbouring three-site blocks 
(Fig. 1). We illustrate the scheme of the RG calculations using the recursion relations 
(3) in the Appendix. 

Using these recursion relations one can find out the ground state energy Eq by 

Eo = hm + NU/2 , 

n— >oo 

where the superscript refers to the n-th stage of iteration. The term NU/2 in the above 
expression comes from the constant term added in (2). One could study the recursion of 
any suitable operator to calculate any other physical quantity within this approximation. 
For example, we calculate the local moment Lq given by 

^0 = ^ (n| - nif 
The recursion relation for Lq turns out to be 

Lo = A + BL'o 

where A and B are numbers depending on the components of the basis states retained 
in course of basis-truncation. Since all the three sites in the basic triangular block are in 
the same status (after adapting the C^y symmetry) one can use any of them for finding 
the recursion of Lq. 



It is interesting to note that due to the group theoretical restriction on the choice of 
the truncated basis it is not possible, in general, to select the lowest energy states from 
all the four groups at a time. The states with the lowest energies in all the four groups 
do not necessarily belong to the same irreducible representation of C^y. If such states are 
retained then the matrix element A defined in (4) will be zero. So instead of the states 
with lowest energy one targets the states belonging to the same irreducible representation. 
This problem was not present in the 1-D cases. Therefore, it is possible here to choose 
any of the possible combinations of four states (from the four groups) compatible with the 
symmetry group. Again from (3) one can see that the contribution to the ground state 
energy from each iteration is oc (E2 + Es). Then it is natural to seek for the combination 
which gives the lowest value of this quantity. However, one should be careful about the 
value of the renormalized hopping t' thus generated, because this in turn will seriously 
affect the contribution to the energy in the subsequent iteration. We have, therefore, 
taken into consideration all the possible channels permitted by the symmetry group up to 
three subsequent iterations and then considered the energetically best ones for the next 
iterations. For example, consider that we have an optimum choice of five distinct channels 
(each with a unique set of four states belonging to a given irreducible representation) to 
start with. In three successive iterations it gives rise to 125 (=5x5x5) channels out of 
which we retain the best twenty five for the next step. This is an optimized way to achieve 
the true ground state. Different possible combinations of the two- and three-particle states 
are shown in Table 1 in terms of the irreducible representation, and corresponding energies 
E2 and E^. It is evident from the table that the contribution [E2 + E^)/2 to the constant 
term D in the Hamiltonian (2) makes only a few combinations preferable in constructing 
the RG scheme. 
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Table 1: Energies of the eigenstates of the block-Hamiltonian corresponding lov — 2 and 
V = 3 are E2 and £"3 respectively. These are tabulated for different possible combinations 
of states belonging to different irreducible representations of C^^. We adopt standard 
group theoretical notation [18] {A^ for 1-D irreducible representations while E for the 2-D 
one). 

III. Exact Results for the Model 

Exact ground state energy could be calculated for the present model (2) in two opposite 
extremal ranges of the parameter space in U/t. To do this we try to match a variational 
upper bound E'up to the exact energy to a lower bound E\o. To find out the upper bound 



Eup we choose a trial state | ^triai) and calculate, by variational principle, 
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Eup — 

To calculate the lower bound £"10 wc break up the Hamiltonian (2) into a sum of the 
cluster-Hamiltonians corresponding to the smallest triangular clusters. These could be 
exactly diagonalized. Then, again by variational principle, Eio is given by 

a 

where -EJ^in ^^e lowest eigenvalue for the a-th cluster and C is the constant term, if any, 
appearing in the Hamiltonian. To find out the lower bound we rewrite the Hamiltonian 
(1) in the following form 



27V 



a=l 



(ij}ea,a iea 



4-^ 

2N jj 

= E ^cluster(a) + -A^ , (5) 
a=l ^ 

where we have used n = U/2 and compensated for the same by adding a term NU/2, 
N being the number of lattice sites. So, here, C = NU/2. In the above form of the 
Hamiltonian we have summed over all possible (2A^ in number) triangular clusters a. 
The fractional numbers appearing with t and U are due to the over counting of the bonds 
(twice each) and the sites (each being shared by six adjacent triangles). i?ciuster(Q!) is the 
cluster Hamiltonian of the a-th cluster for which we have to find out the lowest eigenvalue 

mm 

We categorize the results in two regimes as follows: 
1. U > 0: For very large positive values of U the system is expected to go over to a phase 
with singly occupied sites [17]. We choose the trial wavefunction as 

m = U4,u<^]i io), (6) 

ieC jeC 

where | 0) is a site-vacuum and C and C are arbitrary disjoint sets of lattice sites, each 
containing N/2 lattice points, which together build up the whole lattice. Using this l^f) we 
obtain an upper bound E^^p — 0. We find that the different energy eigenvalues belonging 
to different values of z/ in a basic triangular plaquette a to be as shown in Table 2. Now 
the lowest of these, -E^im "^i^^ be —U/A if 

- C//4 < Min [±t - U/12, ±t/2 - U/12, -t - U/6, t/2 - U/Q, 0] . (7) 

This is true for t/ > 12 | t |. Then E^^ = -U/i. Consequently E^^ = 2iV x (-t//4) + 
NU/2 = = E,p. 

So, in the U > sector exact ground state energy £"0 = for [/ > 12 | t | and the 
ground state becomes a paramagnetic insulator with all sites singly occupied. 
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Table 2: Eigenvalues of the cluster-Hamiltonian i?ciuster(Qi) as given in (5) corresponding 

to different number (z/) of particles. 

11. U < 0: Here we choose the trial wavefunction as follows: 

m = l[4Ai io)' (8) 

where C denotes a set of arbitrarily chosen N/2 sites. This choice gives E^p = NU/2. 
Now using the possible values of the lowest energies in a triangular plaquette as listed in 
Table 2 we find that £^min = if 

< Min [±t - U/12, ±t/2 - U/12, -t - U/6, t/2 - U/6, -U/A] . (9) 

This corresponds io U < —12 \ t \ and consequently Eio = NU/2 = E-^p . 

So, in the U < sector exact ground state energy Eq = NU/2 for U < —12 | t \ and 
the ground state becomes an insulator composed of N/2 "doublon"s. 

IV. Results obtained in the RG 

The ground state energy as a function of the coupling constant U/t, calculated in the 
RG, is shown in Figs. 2(a) and 2(b) for the two cases, t > and t < respectively. It 
appears that the curves completely agree with the asymptotic solutions given above. In 
fact Eq reaches the zero value well before the threshold value of C/, C/ci = 12 | t | as one 
increases the value of U/t. However, there appears a difference between the two values 
of C/ci for t > and t < 0. This merely reflects the lack of t ^ — t symmetry owing to 
the tripartite nature of the lattice. Also on the other side of f//t = the energy curves 
are in complete agreement with the available exact results. Here also Eq takes up the 
value NU/2 at a much higher value of U compared to the threshold value Uc2 = ~12 | t | 
as the value oi U/t is decreased. It appears that the wide region in the parameter space 
(—12 \ t \< U < 12 I i I) for which the exact solution could not be obtained is much 
narrower in reality. The energy curve obtained in the RG calculation smoothly interpolates 
between the two exactly solvable opposite limits. In the intermediate region we find that 
t > always gives the lower energy. This can be checked against a naive calculation for 
any reasonable size of a cluster having a triangular geometry. 

The local moment, which measures the proliferation of "doublons" in the ground state, 
is plotted against U/t for t > and t < in Figs. 3(a) and 3(b) respectively. As we can 
readily see the singly occupied insulating phase given by the wavefunction (6) corresponds 
to Lo = 3/4 = 0.75 lor U > 12 \ t \. Similarly for U < -12 | t |, the wavefunction (8) 
corresponds to Lq = 0. These are reproduced in the RG in accordance with the energy 
curves. The parameter space in between these two insulating phases shows a non-trivial 
discrepancy between the cases t > and t < 0. A large plateau at Lq = 0.5 appears for 
t > which is totally absent in the case of t < 0. For the latter we find a wide plateau 
at Lq = 0.25. Other small structures in the local moment curves may be consequences of 
finite size effects. The plot of the local moment clearly shows that the extent of pairing 



(in the form of doublons) is different in the two cases although both of them possibly 
correspond to a metallic ground state. It is important to note that such a metallic phase 
in the same model on 1-D chain gives free fermionic local moment Lq = 0.375 [6]. These 
plots also show that the phase transitions occurring at 1/^ and are abrupt as in the 
ID case [6]. The departure from Lq = 0.375 in the metallic case is a fall-out of the lattice 
geometry. To naively illustrate this point let us consider the specific case oiU —\ t \— 1.0. 
As one can readily check from Table 2 the lowest energy in each triangular plaquette a 
comes from z/ = 2,4 if t > 0; typical configurations corresponding to this is shown in 
Fig. 4(a). Clearly such configurations will dominate in the global wavefunction and the 
value of Lo will be pushed towards 0.5. On the other hand, a similar observation reveals 
that in case oi t < lowest contribution comes from the i/ — 3 sector with configurations 
similar to those shown in Fig. 4(b). These obviously lower the value of Lq towards 0.25. 
In reality, however, Lq is slightly greater than 0.25 at this point for t < as one can see 
from Fig. 3(b). This is because of the mixing of other configurations for optimization of 
the hopping process between different clusters. Also the finite-size effects might be there. 

The difference in the values of Uc^ for i > and t < is distinctly visible from the 
plots of the local moment. Thus the local moment plot very clearly captures the lack of 
t ^ —t symmetry which is essential for the lattice under consideration. 

V. Conclusion 

Summarizing, we have studied the Hubbard model with the bond-charge interaction 
on a triangular lattice at the special point of particle-hole symmetry. At this point we 
obtain exact results in two opposite limits of the parameter space. The system behaves 
as paramagnetic insulator above a certain value Uc^ of the on-site correlation. Below the 
critical value it undergoes a transition to an insulating state of disordered doublons. 
To explore the full parameter space we employ a real space version of the RG. The RG 
scheme is suitably adapted for this purpose. The ground state energy and the local 
moment values calculated in the RG reproduce the exact results as it did in the 1-D case 
[6] too. This lends some support for the present RG approximations. In the intermediate 
range of the parameter space, where no exact solution has been available, the RG results 
indicate that the degree of pairing (in terms of the formation of local doublons) is different 
from that in the 1-D counterpart of this problem. The triangular geometry plays a crucial 
role there. Moreover, in this region, both the energy and the local moment plots show 
up the effect of loosing the t ^ —t symmetry. The parameter space between Uc^ and 
Uc2, corresponding to a possible metallic phase (as it was in the 1-D counterpart [6, 19]), 
appears to be less wider in the RG calculations compared to the exact solution. The 
present study gives an indication that the RG scheme used here could be successfully 
used to other cases in 2-D. Extension to the cases X ^ t seems interesting although 
it is well known that the lack of particle-hole symmetry creates some problem with the 
present form of the RG on a non-bipartite lattice. It is interesting to look for the short- 
ranged correlations, if any, in the intermediate region of the parameter space. It is also 
interesting to know the effect of the finite block-size on the satellite plateaus in the local 
moment plot; this requires a larger block in the RG analysis (a body-centered hexagon is 
the choice next to a triangle). However, the convergence of the RG results to the exact 
ones is often slow with the increasing block-size, and the effect of the discarded states 
may be important in determining the global wavefunction [8, 20]. Therefore, a better way 
of supplementing the present study is to use the DMRG algorithm, which can take into 



account a larger number of configurations within a block in a controlled and systematic 
way. Of course the DMRG algorithm has to be suitably adapted (on a 2-D lattice in the 
thermodynamic limit) for this purpose. 
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Appendix 

In the present appendix we give a brief sketch of the derivation of the recursion relation 
(3) quoted in the text. We have retained a suitable state of low energy in each of the 
subspaces {u — 2, S — = 0}, {u = 3, S = = |}, {v = 3, S — —S^ — |} and {u — 
A, S — Sz — 0}. These states are designated by 1 0'), |t'), and |tJ,') respectively. Let the 
corresponding energy eigenvalues be Eq/, Ey, E^i, E^i respectively. As we have mentioned 
earlier the first and the fourth states are connected by the particle-hole symmetry while 
the second and the third, by spin-reversal symmetry. It follows, therefore, that Eq/ = 
E-^y — E2 and Ey — Ei' = E3. These four states closely resemble the single-site states 
I 0), It), \l) and It J,) in that the spin quantum numbers are the same and there is a 
one-to-one correspondence between the total number of electrons, u, in such a state and 
the occupation number of the corresponding single-site state (they differ by two). We, 
therefore, identify a three-site block as a "renormalized" site in the scaled lattice and the 
retained states of the three-site block as the "renormalized" single-site states [3, 4]. 

The intrablock Hamiltonian could be written within the subspace of the truncated 
basis in terms of the new block-fermion operators 

H'^ = ^{E, + E,) + 2(^2 - E,) Q - n^>i Q - n^,) (10) 

where the prime denotes the renormalized block-fermion operators. Comparison of this 
"renormalized" block Hamiltonian with the single site part of the Hamiltonian (2) leads 
to the renormalization formulae for U and D in (3). 

To obtain the interblock part of the Hamiltonian we calculate the matrix elements of 
the old fermion operators on the boundary site h, c^i and C(,|n6|, between the states we 
have retained. This leads to the renormalization of t as follows 

t' = 2(AA* - AA'* - A'A*) (11) 

which is equivalent to the last relation in (3). A and A' are already given in (4). 

Let us now illustrate the scheme of calculation of the renormalized parameters as 
given by Eqn. (3). We refer to a specific case with a choice of states belonging to the 
^'-representation (Table 1) from all the subspaces. We find from the Table 1 that there 
could be several sets of states from the subspaces with f — 2 and v — 3 i.e. with 
eigenvalues E2 and £'3 respectively. For example, we pick up the case with E2 — t — U/A 
and ^3 = -3U/A. 

The block-state |0') (with 6*2 = and v = 2) which is an eigenvector of the block- 
Hamiltonian belonging to the eigenvalue E2 = t—U/i and is simultaneously an eigenvector 
of the rotation operators R{C3y) is given by 

|0') =-^[|0 t i)-|0 i t)-2|t i) + 2|i t) 



+ IT i o)-|i t o)] + ^[|o T i)-|o i T) 

-IT i 0) + |o i T) ] . (12) 

Similarly, the block-state | f) (with — 1/2 and u — 3) which is an eigenvector of 
the block- Hamiltonian belonging to the eigenvalue — — 3C//4 and is simultaneously an 
eigenvector of the rotation operators R{C^y) is given by 

IT') =^[lt i T)-2|T T i) + U T T)] 

+^ [ I T i T> - U T T) ] ■ (13) 

Here we have expressed a given configuration of the triangular block by \f3 7 5) {f3, 7, 5 = 
0, T, i or 11) such that the configuration obeys the ordering of the site indices 1 ^ 
2 — >• 7 and 3 — >• 5 (see Fig. 1 for the site indices). 
For the choice of these particular set of states, 

A = (o'|ciTir) = -r^ + 



2V2 2V6 ' 

while 

y=(0'|ci^nu|T') = 0. 

A' may have non-zero value for other sets of states. Instead of using and one could 
use operators belonging to site no. 2 or 3 to land up with the same result. This is because 
all the three sites are "boundary sites" (as pointed out earlier) and, therefore, equivalent 
to each other. These lead to the following RG equations 

U' = 2(^2 - E3) = 2t + f/ 

D' ^3D+ {E2 + ^3)/2 = 3D + t/2 - U/2 

t' = 2Re[ A*(A - 2A') \ t^^t . (14) 

The renormalized parameters could be easily calculated in a similar way for any other 
choice of states shown in Table 1. 
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Fig. 1 Two neighbouring triangular blocks used in BRG. Sites are numbered in a 
counter-clockwise direction. Broken lines show the two connecting paths between 
the two blocks. One has to calculate contributions (identical for symmetry-adapted 
blocks) from both the bonds in calculating the renormalized hopping (3). 
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Fig. 2 Plot of Eq/N, the ground state energy per site for the Hubbard model 
with the bond-charge interaction X = t on a triangular lattice, as a function of 
U/t for (a) t > and (b) t <0. The solid curves show the BRG results while the 
dotted lines correspond to Eo = NU/2 for U/t < and Eo = for U/t > and 
thus show the asymptotic exact results in respective limits. 
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Fig. 3 Plot of Lq, the local moment for the Hubbard model with the bond-charge 
interaction X = t on a triangular lattice, as a function of U/t for (a) i > and 
(b) t < 0. 




Fig. 4 Energetically most favourable configurations in a three-site cluster for 
a typical point in the metallic region, U =\ t \= 1. (a) Configurations with 
u = 2,4 dominate for t > 0. Such configurations correspond to Lq = 0.50. (b) 
Configurations with v = 3 dominate for t < corresponding to Lq = 0.25. In this 
schematic figure f (J,) stands for an electron with "up spin" ("down spin") at a 
given Wannier orbital, o denotes an empty site. 



